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There is much current interest in modelling suspensions of algae and other micro- 
organisms for biotechnological exploitation, and many bioreactors are of tubular 
design. Using generalized Taylor dispersion theory, we develop a population-level 
swimming-advection-diffusion model for suspensions of micro-organisms in a vertical 
pipe flow. In particular, a combination of gravitational and viscous torques acting on 
individual cells can affect their swimming behaviour, which is termed gyrotaxis. This 
typically leads to local cell drift and diffusion in a suspension of cells. In a flow in 
a pipe, small amounts of radial drift across streamlines can have a major impact on 
the effective axial drift and diffusion of the cells. We present a Galerkin method to 
calculate the local mean swimming velocity and diffusion tensor based on local shear 
for arbitrary flow rates. This method is validated with asymptotic results obtained 
in the limits of weak and strong shear. We solve the resultant swimming-advection- 
diffusion equation using numerical methods for the case of imposed Poiseuille flow 
and investigate how the flow modifies the dispersion of active swimmers from that 
of passive scalars. We establish that generalized Taylor dispersion theory predicts 
an enhancement of gyrotactic focussing in pipe flow with increasing shear strength, 
in contrast to earlier models. We also show that biased swimming cells may behave 
very differently to passive tracers, drifting axially at up to twice the rate and diffusing 
much less. 
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I. INTRODUCTION 

Swimming micro-organisms, such as algae and bacteria, have their own agenda; selec- 
tive pressures lead cells to adopt strategies to optimize a combination of environmental 
conditions, such as illumination, nutrients or the exchange of genetic material. This can sig- 
nificantly impact the behaviour of suspensions of swimming micro-organisms, particularly 
in flows where biased motion across streamlines can lead to rapid transport. For example, 
various algae are gravitactic, that is they swim upwards on average in still fluid which can 
be beneficial for reaching regions of optimal light. For some species this is due to being 
bottom- heavy - the centre of gravity for these cells is offset from the centre of buoyancy, and 
the combination of the effects of gravity with the buoyancy force gives rise to a gravitational 
torque which serves to reorient the cell allowing it to swim upwards - whereas in others 
sedimentary torques lead to similar behavioui'^. However, in shear flow the cells may be re- 
oriented from the vertical due to viscous torque^. For a vertical pipe containing downwelling 
fluid, gravitactic cells can accumulate near the centr^, a phenomenon known as gyrotactic 
focussing. As recently predicted theoretically by Bees and Croze^l, such a modification of 
the spatial distribution of algae in tubes alters significantly the effective axial dispersion of 
the cells. 

There is much current interest in employing micro-organisms for biotechnological pur- 
poses, from the production of biofuel^^, such as hydrogen, biomass or lipids, to high-value 
products, such as /3-carotene. Cells are grown either extensively on low value land or in- 
tensively to optimize growth. Intensive culture systems typically consist of arrays of tubes 
(vertical, horizontal or helical) and aim to maximize light and nutrient uptake. Bioreactors 
may be pumped or bubbled, in turbulent or laminar regimes. However, energy input may be 
energy wasted; efficient bioreactor designs might aim to make use of the swimming motion 
of the cells themselves, or accommodate the fact that swimming micro-organisms (where 
drift across streamlines is more important than axial motion) and nutrients are likely to 
drift and diffuse at different rates along the tubes. 

In a still fluid, the swimming behaviour of individual gyrotactic phytoplankton has been 
usefully described as a biased random walk: the cell orientation is assumed to be a ran- 
dom variable that undergoes diffusion with drift^. At the population-level the dynamics can 
be modelled with a swimming-diffusion equation for the cell concentration, where the cells 
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swim in a preferred direction at a mean velocity and diffuse with an anisotropic diffusion ten- 
sor that represents the random component of swimming. Extending such population-level 
models to incorporate the effects of ambient flow is non-trivial. Although the orientation 
distribution and resultant mean swimming velocity of such cells in unbounded homogeneous 
shear flow has previously been computecP'^, the resultant diffusion tensor is more compli- 
cated. For homogeneous shear flow, subject to certain constraints on the form of the flow, 
Hill and Bees'^ and Manela and Frankell^ calculated expressions for the diffusion tensor 
using the theory of generalized Taylor dispersion (GTD). Because of its account of shear- 
induced correlations in cell position, GTD is a more rational account over earlier approaches 
based on an orientation only description using a Fokker-Planck equation and diffusion tensor 
estimate (FPf^. 

Bearon, Hazel, and Thornf^l compared two-dimensional individual-based simulations of 
swimming micro-organisms with swimming-advection-diffusion models for the whole popu- 
lation in situations where the flow is not homogeneous, that is in flows in which the cells can 
experience a range of shear environments. Using GTD theory to calculate local expressions 
for the mean swimming direction and diffusion coefficients, the results of the individual and 
population models were generally in good agreement and were able to successfully predict 
the phenomena of gyrotactic focussing. However, this work was restricted to two-dimensions; 
both the swimming motions and velocity field were confined to a vertical plane. 

Here, we consider axisymmetric pipe flow, which locally can be described by planar shear, 
and consider swimming motions which are allowed to be fully three-dimensional. First, we 
develop a population-level swimming-advection-diffusion model where the mean swimming 
velocity and diffusion tensor are based on the local shear. Next, a Galerkin method is 
presented for calculating the mean swimming velocity and diffusion tensor based on the 
local shear, and asymptotic results are obtained in the limits of weak and strong shear. The 
resultant swimming-advection-diffusion equation is then solved numerically for the case of 
imposed Poiseuille flow. We contrast the GTD results with the FP approach. Finally, we 
investigate how the flow modifies qualitatively and quantitatively the dispersion of active 
swimmers from that of a passive scalar. 

This paper represents an important link study that will facilitate the comparison of the 
exact long-time theoretical results of Bees and Croze'^ and the forthcoming experimental 
results by the authors on the transient dynamics. 
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II. MATHEMATICAL MODEL 

A. Vertical pipe flow 

Consider axisymmetric fluid flow with velocity u through a vertical tube of circular cross- 
section, radius a, with axis parallel to the 2;-axis pointing in the downwards direction, such 
that 

u = u{r)e^ = U{1 + xir/a))e^. (1) 

Here, U is the mean flow speed, Ux is the variation of the flow speed relative to the mean, r is 
the radial distance from the centre of the tube and (e^, e^, e^) are right-handed orthonormal 
unit vectors that define the cylindrical co-ordinates. For flow subject to a uniform pressure 
gradient and no-slip boundary conditions on the walls, we have simple Poiseuille flow, x(r) = 
1 — 2r^. In the fully coupled problem, where the negative buoyancy of the cells modifies the 
flow, x(r) must be determined, as in Bees and Croze'^. 

A population-level model for gyrotactic micro-organisms in homogeneous shear flow has 
previously been derived based on generalized Taylor dispersion theorj^^^ (GTD). Specifl- 
cally, for particular types of flow and on timescales long compared to 1/dr, where dr is the 
rotational diffusivity due to the intrinsic randomness in cell swimming, the cell concentration 
n{'x,t) was shown to satisfy a swimming-advection-diffusion equation of the form 



+ Vx. 
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where Vg is the constant cell swimming speed, and q and D are the non-dimensional mean 
cell swimming direction and diffusion tensor, respectively. Explicit expressions for q and 



D as a function of the local shear strength will be given in section II B Furthermore, 
Bearon, Hazel, and ThomfHl show that this population-level approach is a good approxi- 
mation for flow flelds more general than homogeneous shear. Therefore, we shall use ^ 
to describe the cell concentration in a pipe flow with non-homogeneous shear. To solve the 
swimming-advection-diffusion equation numerically, it is convenient to non-dimensionalize 
lengths based on the pipe radius, a, and non-dimensionalize time on o?dr/V^, a characteris- 
tic timescale for diffusion across the pipe. This reveals two non-dimensional parameters in 
the problem: the Peclet number which is given by 

Pe = ^, (3) 
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and /5, the ratio of pipe radius to a typical correlation length-scale of the random walk in 
the absence of bias, defined as 

/'-f. (4) 

An alternative interpretation of /3 = aVs/ {Vgd'^) is as a 'swimming Peclet number'. Equa- 
tion ([2]) in non-dimensional form thus becomes 

— + Vx. [{Pe[l + x(r)]e, + /3q)n - Ti.V^n] = 0. (5) 

B. Generalized Taylor dispersion 

The shear in the pipe flow given by ([T]) can be locally described as a simple shear flow. 
Specifically, consider a Taylor expansion of the flow field near some reference point Rq which 
is at radial position r = Rq 

u(R) ^ u(Ro) + (R - Ro).e,-x'(i?o/a)e,. (6) 

CL 

We consider local co-ordinates relative to an origin located at Rq such that k is pointing 
vertically upwards and (i, j, k) form a right-handed orthonormal set of unit vectors so that 

i = e^, j = -e^, k = -e^. (7) 

Defining the local position co-ordinate, R — Rq = ■Ci + '7j + Ck, the flow field can then be 
written locally as simple shear, such that 

u(R) = u(Ro) + G^k, (8) 

where the shear strength G is given by — ^x'- With this choice of co-ordinates, the velocity 
gradient tensor, G, defined such that u(R) = u(Ro) + (R — Ro).G, has the simple form 
Gij = G6ii6j3. 

The mean swimming direction, q, and non-dimensional diffusion tensor D can be written 
as integrals over cell orientation, p, in the formfnUSl 

q = / Pf{p)dp, (9) 



p 



D = / [bp + -^bb.G]"^"t;p. (10) 



p 
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Here [j'^^"^ denotes the symmetric part of the tensor, G = ik, and o" is a non-dimensional 
measure of the shear, defined as 

G Pe , 

" = 2d; = -2^^- 

We note that a varies with r because the shear varies across the radius of the tube. However, 
in the theory of GTD the shear is assumed locally homogeneous, and so we calculate local 
expressions for the mean swimming and diffusion based on the local value of a. 
The equilibrium orientation, /(p), and vector b(p) satisfy*^^^ 

Cf = 0, (12) 
£b-2ab.G = /(p)(p-q), (13) 

subject to the integral constraints 

fdp = l, I hdp = 0. (14) 
Here, the linear operator C for a spherical swimming cell is defined by 

Cf = Vp.((A(k - (k.p)p) - aj A p)/ - Vp/), (15) 
the gyrotactic bias in swimming direction is represented by the non-dimensional parameter 

and B = fia±/2hpg is the gyrotactic reorientation time scale, where h is the distance between 
an average cell's centre-of-mass and centre-of-buoyancy, a± is the dimensionless resistance 
coefficient for rotation about an axis perpendicular to p, /i and p are the fluid viscosity and 
density respectively, and g is the magnitude of the gravitational force. 

To summarize, the non-dimensional mean swimming velocity and diffusion tensor are 
given as functions of two non-dimensional parameters: A, which only depends on proper- 
ties of the cell, and a, which quantifies the strength of the shear. See equations (|9 16). 



Furthermore, for the pipe flow considered in the previous section, we can express a as a 
simple function of the non-dimensional parameters Pe, the global Peclet number, and /3, 



the swimming Peclet number, and the shear profile x'(^) (Eq. 11). The solution of the 
governing non-dimensional swimming-advection-diffusion equation (Eq. |5]) can therefore 
be determined by specifying the three non-dimensional parameters A, Pe and /3 and the 
non-dimensional flow proflle x(^)- 
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When explicit calculations are presented in this paper we have assumed that the flow is 
Poiseuille, x{i^) = 1 ~ 2r^, and take A = 2.2 so as to compare with previous worliP^ based 
on the algal species C. augustae (wrongly identifled as C. nivali^). In Bearon, Hazel, and 
Thornl^, good agreement was found in planar pipe flow between individual based simulations 
and the population- level model with /3 = 10 (where the reciprocal of /3 was deflned as e = 0.1 
therein). Motivated by the pipe dimensions in experiments currently in progress, we also 
consider (5 = 2.34. Note that /3 represents the ratio of pipe radius to the correlation length- 
scale of the random walk in the absence of bias. Therefore, when modelling a random walk 
as a diffusion process, /3 should be sufficiently larg^. However, we hypothesize that this 
restriction may be relaxed in the case of gyrotactic cells that are well-focussed by the flow 
along the axis of the tube and only suffer rare collisions with the wall. 



C. Calculation of mean swimming velocity and diffusion 

Hill and Bees'^ demonstrate that the GTD equations (12) and (13) for / and b, re- 



spectively, can in general be solved by expanding in spherical harmonics using a Galerkin 
method. The method is summarized for the flow employed in this paper in appendix |Xj 

To simplify the numerical solution of the swimming-advection-diffusion equation in pipe 
flow, we flt the rather complex algebraic expressions in a obtained using the Galerkin method 
for the mean swimming direction and diffusion tensor with the simpler curves 

g'^(a) = -aP(a;a^bO, g^(a) = -P(a; a^ b^), (17) 
D^ ia) = P{a- a^^ b^), D'-{a) = -aP{a- a^, b^), D^^{a) = P{a; a^, b^), (18) 

where 

The choice of a and b coefficients is described in table [T] with reference to asymptotic results 
presented below. Please refer to appendix |E] for the coefficients of fits to the full Galerkin 
solution. We also consider results using the simpler estimate for the diffusion tensor, which 
we describe as the Fokker-Planck approximation (or FP), discussed in detail in appendix [Pj 
In flgure [T] we see that these simple functions are good approximations for the exact 
solutions for the mean swimming and diffusion, and it is evident how shear can signiflcantly 
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TABLE I. In order to obtain the simplest functional fits whilst ensuring the asymptotic results 
are satisfied, the a coefficients are as specified. The free parameters, 62)^4 are obtained through 
least squared optimization of each fit of velocity or diffusion component against a. Because we are 
unable to obtain easily the coefficient of the 0{a) correction to D'q , in addition we allow a^^Q to 
vary. The fit coefficients for A = 2.2 are given explicitly in appendix [E] The subscript G highlights 
that the results are for generalized Taylor dispersion (GTD), and ** indicates that the parameter 
is fitted. 
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affect the mean swimming and diffusion. Furthermore, we note how that the calculation 
of diffusion via the GTD method is qualitatively different to that calculated via the sim- 
pler Fokker-Planck method (FP). In particular, we note that the components of diffusion 
approach zero in the limit of large shear using the GTD method, whereas they approach a 
finite non-zero limit via the FP method (see Hill and Bees^. 

To provide confidence in the results from the Galerkin method, we have obtained asymp- 
totic expressions for cr ^ 1 and a ^ 1, as described in appendices [B] and [C| 

Specifically, for cr ^ 1, the mean swimming direction with respect to coordinates 
(e^, e^, e^) correct to 0((t) is given by 

q=-{jJ,,0,K,f, (20) 

where the quantities Ji and Ki are specified functions of A, coinciding with the results 
of Pedley and Kessler^^ using the FP model. However, calculation of the diffusion tensor 



from (10) reveals the new result that at leading order the diffusion tensor is diagonal with 



horizontal component D"' = jh, and vertical component D^^ = where Li is also a 



specified power series in A (appendix B2). For A = 2.2 we have that Ki = 0.57, Ji = 0.45 



and Li = 0.11, see appendix|B| There is an 0(a) correction to the off-diagonal term D^'^, but 
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FIG. 1. Mean swimming and diffusion coefficients as a function of shear, a, for A = 2.2. Points 
are calculated using the Galerkin method, solid lines are functional fits described in the text, and 
dashed lines are asymptotic results. For diffusion calculations, black lines are for GTD, whereas 
red (grey) lines indicate the FP estimate. 



the second term in the definition of the diffusion tensor in ( 10 ) does not allow for a simple 
closed form expression for these components. (This is in contrast to expressions obtained by 
Pedley and Kessler'^ using the simpler orientation-only FP model with a diffusion estimate 
proportional to the variance of p.) 

For 0" ^ 1, as in Bees, Hill, and Pedleyl^, the mean swimming direction correct to 0(l/cr^) 

is 

/ 2A 4A .rr, / ^ S 

Here, using GTD theory, we have the new result that the non-zero coefficients of the diffusion 
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tensor are 



D 



v 



„2 



d2 
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3^ ^ / 



4 + ^ 



(22) 



where the quantities di, c?2, d^, d^, d^ are polynomials in A, given in appendix [C| For A = 2.2, 
we find that d^ = 0.68, di = 0.0060, d^ = 0.0020, d^ = 5.9, 4 = 1.3. 

The asymptotic results are presented in figure [T| indicating excellent agreement with 
results from the Galerkin method and demonstrating correspondence with the functional 
fits described above. 



III. POPULATION-LEVEL NUMERICAL SIMULATIONS 
A. Numerical Methods 



The governing swimming-advection-diffusion equation is solved using a spatially adaptive 
finite element method as described in Bearon, Hazel, and Thorn The cell concentration, n, 
is approximated using standard Lagrangian quadratic finite elements and the time derivative 
is approximated using an implicit second-order, backward difference scheme. The subsequent 
discrete linear system is assembled using the C++ library oomph-lib'^and solved by a direct 
solver, SuperLU^^. In unsteady simulations, a fixed time-step of dt = 10~^ is used. The 
results were validated by repeating selected simulations with smaller error tolerances and 
timesteps. 



B. Steady gyrotactic focussing 

First, we seek an equilibrium solution n(r) of equation ^ which represents gyrotactic 
focussing of cells towards the centre of the pipe. Imposing zero fiux on the pipe wall, at 
r = 1, we have that 

dii 

/Sq^n - D'^^ = 0, (23) 
dr 

which we can integrate to obtain 

n = noexp(^j j^dr^ , (24) 
10 
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where the radial components of the mean swimming direction and diffusion tensor, and 
D^*", respectively, are functions of the local shear. In particular, if we take simple Poiseuille 
flow, x(r) = 1 — 2r^, we have that a, the non-dimensional measure of the shear, is given by 
a = —^x' = yrr. For a ^ 1, at leading order we have that q'' jD" = —aX = —^^r 
from which we predict the Gaussian distribution 

n = noexpl —r \. (25) 

As demonstrated in figure [l](f), the leading order asymptotic solution q^ / D^^ = —a\ is 
an excellent approximation for a = 0(1). It is important to note that the GTD and FP 
methods yield a qualitative difference in the behaviour of q^/D''^. 



Example calculations of the equilibrium solution (Eq. 24 ) are shown in figure |2} For 
the given values of Pe and (3, we see that cells undergo gyrotactic focussing, and that the 
distribution predicted by GTD theory can be well-approximated by the Gaussian distribution 



(Eq. 25) but shows a marked difference to that predicted by FP theory. Furthermore, 
whereas we see that GTD predicts an enhancement of gyrotactic focussing with increasing 
shear strength, the FP approximation predicts a reduction in gyrotactic focussing with 
increasing shear strength at sufficiently large shear. 



C. Vertical dispersion 

Bees and Croze"^ investigated how the average axial dispersion was modified for gyrotac- 
tic organisms compared with a passive solute. Specifically, using the method of moments 
and the FP approach they obtained long-time expressions for the vertical drift relative to 
the mean flow and the effective axial swimming diffusivity as a function of Pe and a gyro- 
tactic parameter. Here, we perform a similar calculation, using simulations and the GTD 
calculations for the diffusion tensor. We solve numerically the swimming-advection-diffusion 
equation ([s]) with initial condition 

«(r...O) = n„exp(-(i^)'-(i^)'), (26) 

representing a Gaussian blob of cells centred at z = O.IL, r = 0. For the simulation domain 
we take 2; G (0,L), r G (0,1). Furthermore, we impose no-flux boundary conditions on 
the walls r = 1, symmetry around the centreline, and periodic boundary conditions in the 
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FIG. 2. Equilibrium concentration (24) for Pe = 20 (solid line) and Pe = 50 (dashed line) for 
swimming parameter /? = 2.34. Cell diffusion is calculated using GTD (black) and FP (red/grey) 



approaches. The dotted lines are the associated Gaussian distributions (25). The solutions are 
normalized so that there is unit total mass per unit length, J 2TTrn{r)dr = 1. 



vertical direction, but take L to be sufficiently large that boundary effects do not influence 
the vertical distribution. In the results presented in figures [3] and |4| we take Pe = 50, 
L = 1200 and run the simulations for t G [0,8]. In figure |3] we see example plots of the 
early concentration distribution as a function of time for both gyrotactic cells and a passive 
solute. For the passive solute, we take D = |l, and q = in equation (jsj). As shown in 
appendix |B| this is equivalent to considering mean swimming and diffusion in the absence of 
gyrotactic bias, A = 0, and shear, a = 0. Whereas the gyrotactic cells are focussed towards 
the centre of the pipe, the passive solute diffuses radially. 

Following Bees and Croze'^, we quantify dispersion in terms of cross-sectionally averaged 
axial moments of the concentration distribution. To compute the moments of the distribu- 
tion, we ffist translate to a reference frame moving with the mean flow, z = z — Pe t. The 
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FIG. 3. Concentration in region z £ (0,600), r £ (0, 1) from t = to t = 1 at intervals of 6t = 0.1 
with Pe = 50. Upper plots are for gyrotactic cells with A = 2.2, /3 = 10. Lower plots are equivalent 
results for a passive solute. The colour scale is based on the initial concentration distribution, with 
red representing the maximal initial concentration at the centre of the blob of cells, and blue zero 
concentration 

cross-sectional average, mp(t), of the pth axial moment, Cp, is (dropping hats for clarity) 

Cp{r,t) = j zPn{r,z,t)dz, p = 0,1,2, (27) 
mp{t) = 2 j Cp(r, t)rdr, p = 0,l, 2. (28) 

The mean and variance, mi and m2 — vn\, of the distribution are plotted in figure [4| The 
solution is normalised so that the total mass is unity, mo = 1. 

From the calculations of the m\ and m2, we then define the axial drift and effective axial 
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FIG. 4. The mean and variance, mi and m2 — ml of the distribution as a function of time, t for 
Pe = 50. Upper plots are for gyrotactic cells with A = 2.2, f3 = 10. Lower plots are equivalent 
results for a passive solute. Open circles are results from numerical simulation, solid lines are linear 
regressions for t G [4,8]. 

diffusion to be 

Ao = lim ^mi, (29) 

t-s>oo dt 

De = lim ^^(^2 - mj). (30) 
t^oo 2 at 

As depicted in figure |4], for Pe = 50, performing a linear regression over the interval 
t G [4, 8] we obtain Aq = 35.2 for the gyrotactic cells with parameters A = 2.2, /3 = 10, 
compared to the long-time limit of Aq = for a passive scalar predicted from classic Taylor 
dispersion theory. This occurs because gyrotactic cells are focussed towards the centre of 
the tube where the flow is fastest and, hence, they are transported more rapidly than the 
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mean flow. Noting that D = |l for the passive solute, for Pe = 50 the classical Taylor 
dispersion result predicts that Df. = 1/6 + 6Pe^/48 = 313, which compares well with the 
numerical calculation of Df. = 312. For the gyrotactic cells with parameters A = 2.2, 
(3 = 10, we see a much reduced axial dispersion, with an estimate of = 20.0. As 
discussed by Bees and Croze'^, this reduction in axial dispersion can be explained due to 
gyrotactic focussing: by self-concentrating towards the axis of the tube, cells undergo a much 
reduced sampling of radial space and thus sidestep classical shear-induced Taylor dispersion. 
Furthermore, preliminary calculations based on equations 6.1 and 6.2 of Bees and Crozel^ 
using the GTD values for the components of q and D give excellent agreement with these 
numerical computations. Specifically, the calculations yield Aq = 35.2 and Dg = 20.6 for 
gyrotactic cells with parameters A = 2.2, (3 = 10. 

IV. DISCUSSION 

Here, we have considered the spatial distribution of gyrotactic algae in axisymmetric pipe 
flow. We have computed a population-level swimming-advection-diffusion model where the 
mean swimming velocity and diffusion tensor are based on the local shear using the theory 
of generalized Taylor dispersion. We have shown how shear modifies the mean swimming 
velocity and diffusion tensor and, furthermore, demonstrated how the diffusion tensor differs 
qualitatively from previous simpler models, such as the "Fokker-Plank" approach for which 
the diffusion tensor is estimated to be the product of the variance of the orientation distribu- 
tion and a correlation timescale. We have demonstrated that the shear-induced modification 
to mean swimming velocity and diffusion results in gyrotactic focussing and have quantified 
how the axial drift and diffusion of a population of cells is modified from that predicted for 
a passive scalar. 

In this paper, we only considered unidirectional coupling between flow field and cell 
concentration. However, actively swimming cells, that are typically denser than the fluid, will 
modify the flow field. In a dilute suspension, where direct cell-cell hydrodynamic coupling 
can be neglected, negatively buoyant cells modify the flow field from Poiseuille fioA^, which 
results in a change in local shear and thus a modification of the mean swimming velocity and 
diffusion tensor. Furthermore, direct hydrodynamic interactions between cells, and stresses 
induced by the swimming motions, may also alter the flow flelcP^. 
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Work in progress by the authors aims to incorporate the population-level model derived 
here from generalized Taylor dispersion in Bees and Croze'd^l modification of the classical 
Taylor- Aris theory in order to predict the axial drift and diffusion. Furthermore, both these 
predictions of long-time dispersion and the transient results presented in this paper will be 
compared with experimental observations of axial drift and diffusion in dyed suspensions 
of the alga Dunaliella salina in vertical tubes subject to imposed flow. Finally, work is in 
progress by the authors to use direct numerical simulations to study the dispersion of active 
swimmers in laminar and turbulent flows, comparing statistical measures of dispersion from 
simulations with analytical predictions using the GTD expressions derived in this paper. 
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Appendix A: Galerkin method 



To implement the Galerkin method, we follow the approach of Hill and Bees'^ who 



considered the flow fleld u = G^i. Here, for the flow fleld u = C^k (see Sec. HA) we further 



extend the method to establish results for the full positive-deflnite diffusion tensor in (10). 
We parameterize cell orientation in terms of spherical-polar co-ordinates {9, (p) 



p = sin 6 cos (f)! + sin 6 sin (f)j + cos ^k. 



Note that the direction 6 = corresponds to cells directed vertically upwards. 



Equations (12) and (13) are solved by expanding / and bj, j = 1,2,3, in spherical 
harmonics: 



/ = ^^A:rcosm0Pr(cos^), 

ra=0 m=0 
oo n 



(Al) 
(A2) 



n=0 m=0 
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Defining 



F:r,^RnM)Pn{cOSe), J =0,1, 2, 3, 



cos(m0), j = 

/3™ cos(m0) + 7™. sin(m0), j = 1, 2, 3, 



(A3) 
(A4) 



equations ( 12 ) and ( 13 ) then yield 



J2J2{n{n+l)F:[; + X sin^ ^i^^Ti^f' + ^^(cos sin ^i?™ P^' + cot 9 sin 0/?:^/^^ ) 



n=0 m=0 



-2Acos^F„"}} 



Er=oE::.=o(sm^cos0- (47r/3)Al)F-, 

Er=o E^=o sin^sin0F^ 

Er=o E:;=o {2^^n"J + (cos^ - (4vr/3)A?)Fro} 



(A5) 



2, 
3, 



where primes denote differentiation with respect to the dependent variable. Note that the 



normalization condition (Eq. 14) requires that Aq = l/An,(3Qj = 0, and from equation (joj), 
we calculate the mean swimming direction to be q = {4:7!'/3){A\,0, A^)'^ . These equations 
can be simplified using identities for spherical harmonics so that inner products with other 
harmonics can be calculated. Finally, the resulting equations can be approximated by trun- 
cating the above series solutions to give a set of simultaneous equations that may be solved 
for the coefficients A™, /3™ and 7™. 



The first term for the positive-definite diffusion tensor in equation (10), is given in part by 



equation (52) of Hill and Beesl^, which only depends on the first few terms in the expansion 
(i.e. 7]j, for m = 0, 1). The second term cannot be written in such simple terms but can 
be approximated directly using all available coefficients. 



Appendix B: Small a asymptotics 

The calculation of the mean swimming direction, q, is the same for both the generalized 
Taylor dispersion theory and the Fokker-Planck approaclJ^. Hence, we follow Pedley and 
Kessler'^ to compute / for small vorticity case (note that their small parameter e is related 
to a via a = eX). 
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1. Calculation of equilibrium distribution, /, and mean swimming, q 



At leading order, a = 0, (12) becomes 



sm 9 o9 \ oO J sm 6 ocp^ sm 6 oO ^ 
Looking for a solution independent of 0, we obtain the von Mises distribution, 



where (14) yields the normalization constant /i = A/47rsinhA. 



From equation ([9]), the mean swimming velocity is computed to be 



27r /"Tr 

/ p/(°)(^)sin^ci^ci0 = (0,0,iri), 
'0 Jo 

where Ki = cothA — 1/A. For A = 2.2 we have Ki = 0.57. 
To find the 0{a) correction, put 

/ = /(°)(^^) + a/W, 

to obtain 

-^0^^'^ = —n^n 6^^^^) + (^i^'^/^'^) = A cos sin 0, 

smfc'afc' \ o9 J sm 6 acp^ sm9d9 ^ ' 

Looking for a solution of the form f^^^ = cos (f)F(9), we obtain the ODE 
1 d f . dF\ F \ d 



\ d9 J sm 9 sm 9 d9 ^ ' 



sin 9 d9 

Defining x = cos 9, and letting 

we obtain equation (3.4) of Pedley and Kessler^^, which has a power series solution 

OO n 

9li^) = A"^n(a;), An{x) = an,rPr{x), 

n=l r=l 

where the are the associated Legendre functions and the coeffiicents an,r satisfy 

__ r + 2 (r - 1) en+i,r 

(r + l)(2r + 3) r(2r — 1) r(r + 1) 

where 

2r + 1 



e 



n+l,r 



n\2r{r + 1) J_i 
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The first order correction to the mean swimming is then given by 



q 



(1) 



27r /'TT 



^0 



p/^^-* sin 



2lT /'TT 



^0 



p COS (t)F{9) sin OdOdip 



Ji 



^,0,0), 



where 



1,1- 



1=0 



With A = 2.2 a calculation using Maple provides Ji = 0.45. 

Thus the mean swimming correct to 0(cr) with respect to i,j,k unit vectors is given by 
q = {—jJi, 0, . Recalling the relationship between local and global co-ordinate vectors, 



— e 



k = — e. 



the mean swimming direction, with respect to 6^,6^,6^ is (see Eq. 20) 



(B4) 



(B5) 



2. Calculation of b, and diffusion tensor, D. 



At leading order, setting cr = in equation (|13|) yields 
1 d 



-Coh = ^^j^ ( sin^^ 
sm 9 o9 \ o9 



By inspection, consider 



1 d^h \ d 
sin^ 9 sin 9 d9 

= Bh{9) COS0 
6.^ = Bh{9) sin</) 
k = Bv{0). 



(sin^^b) = /(°)(;^ik-p). 



Bh then satisfies the ODE 

1 d . ^ 

sm 9- 



sin 9 d9 



dB 



H 



d9 



B 



H 



A d 



sin 9 sin 9 d9 



{sin^ 9Bh) = -sm9f 



(0) 



(B6) 



On comparing equations (B6) and (B2), we can write 

1 



B 



H 



A 



-F. 



From equation (jlOj), the leading order expression for the non-dimensional horizontal com- 
ponent of diffusion [D^^ = D^^) can thus be written as 



277 /•TT 



JO 



COS (I)Bh{9) sin 9d9d(f) = — 



A 



277 /*7r 



^0 



COS (j)F{9) sin 9d9d(j) = ^ . (B7) 

A 
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The function By satisfies the ODE 

1 d / „dBv\ A d 



sin 9 dO 

To solve this, as for F, we define x = cos 9, let 



and seek power series solutions 



By = jhi{x), 



/ii(x) = ^A"5„(x), 

n=l 
n 



X). 

r=l 

By utilizing properties of Legendre polynomials we obtain the following recurrence rela- 
tionship for the bn,r- 

J bn,r+l ^n,r— 1 fn+l,r 

On+i,r = + ^731 + r(r + l)' 

where 



fn+l,r = J {X- K^)x^P';{x)dx. 

The leading order expression for the non-dimensional vertical component of diffusion can 
thus be written as 

D^^ = i I p^Bv{9)smed9d(l)= (B8) 
Jo Jo A 

where 

4 oo 

Li = — /i^A"6„,i. 

n=l 

For A = 2.2 a computation employing Maple reveals that Li = 0.11. 
The off-diagonal terms, D^'^, etc., are all zero at leading order. 

When comparing with a passive solute, we note that if we set a = 0, at leading order in 
A we have that 

= l/4n, Ji = ^A^ai,i, Li = ^A6i,i. 

Noting that ai^i = 6i i = 1/2, it is clear that in the limit of A — )■ the diffusion tensor tends 
to the isotropic tensor 1/6. This result can be obtained directly by considering equations 



(12 15), which have solution / = l/47r,b = p/Svr. 
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Appendix C: Large a asymptotics 

For the calculation in the limit of large a, it is convenient to follow Manela and Frankel'^ 
and Brenner and Weissman'^ and define local co-ordinates so that the vorticity vector is in 
the direction of k: 

i = e^, j = -e^, k = e^. (CI) 

Defining the local position co-ordinate, R — Rq = + r/j + (^k, the flow field can then be 
written locally as simple shear: 

u(R) = u(Ro) + Gd, (C2) 

where the shear strength G is given as before by — ^x'- Fo^' this flow field, the velocity 
gradient tensor has the simple form Gij = G6ii6j2- 
Writing the orientation vector as 

p = sin 9 cos 01 + sin 9 sin 0j + cos 9k 



we can write the governing equation (12) as 



Cf = a^-Csf = (C3) 
where the linear operator independent of a is given by 

and where in sphericals the Laplacian is given by 

^2r I d f . df\ , 1 d'f 



sin 9d9 \ d9 J sin^ 9 d^"^ ' 
For (T ^ 1, we consider the following perturbation expansions for / and b: 

b=±fbW+-b«+f-Vb(2)+.. 

An \ a V cr / 
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1. Calculation of equilibrium distribution, /, and mean swimming, q 



Substituting the expansion for / into equation ( C3 ) we obtain an explicit iterative scheme 
for computing the expansion: 



d(j) 

At leading order: 

5/(0) 



Csf^'^. (C4) 



0. 



(90 

Hence /(°) = /(°)(^) subject to £ /W(^) sin 6(16 = 2. 

At oc-y. 



^ = CJ<'\6). 



Because f^^^ must be periodic in with period 27r, integrating this equation with respect to 
(j) from to 2tt gives: 

^6d6[''''^^^)='- 
Excluding singular solutions, and given that f^'^\6) sin 6 d6 = 2, we obtain the solution 

We can summarize the general iteration algorithm for the terms k > 1: 



1. Integrate equation (C4) 

/('^+i)= r CJ^''^d<P + F^'^'\6). 
Jo 

2. Impose periodicity of & integral constraint 

Jo 

f^^+^U6d(f) = 0, 



^0 



to determine non-singular solutions for F^'''^^\6). 
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Specifically the first two terms in the expansions are given: 

= -2Acos0sin^, 

/(2) = ^a2(1 - 3 cos^ 9) + 4A sin ^ sin + - cos 20 sin^ 9. 
3 2 

From equation ([9]), we thus can compute mean swimming at large a correct to 0(l/cr^): 

p27V pit r\ \ 

q^= / f sin^ 9 cos (I)d9d(f) = (C5) 
i(/,=o Je=o 3(T 

/•27r A\ 

q"^ = / /sin2^sin0rf^rf0= — , (C6) 

J(f)=Q J 0=0 3(7 

7C 



= / / /cos^sin^c/^rf0 = 0. (C7) 
When converting back to the global co-ordinates we note that 

i = e^, j = -e^, k = e^. (C8) 
and so with respect to e,., 6^,62 unit vectors the mean swimming is given by 

.^-(|.O.i^f. (C9) 

2. Calculation of b, and diffusion tensor, D. 

We now apply similar techniques to calculate the vector field b. Substituting the expan- 



sion for b into equation (13) we obtain the following iterative scheme for computing the 
expansion: 

_ 2b(^+i).G = (47r(p - q)/)(') + Csh^'\ (CIO) 
For the simple shear flow with Gij = SiiSj2, taking the dot product with i yields: 

(47r(sin0cos0-g«)/)^'U£.6f^ (Cll) 

The method follows as for /: 



1. Integrate equation (Cll) 



^(fc+i) ^ (4vr(sin^cos0 - + C,bfd(P + fif 

Jo 
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2. Impose periodicity of 6^'^^^'' & integral constraint 

/ (47r(sin ^ cos - q^)f) ^'^'^ + Csbf-^^dcj) = 0, 
Jo 

Jo Jo 

to determine non-singular solutions for B^^^^\6). 
We obtain the following expressions: 

= 0, 

6 W = A ( 1 - 3 cos^ ^) /36 + sin ^ sin 
bf = Bf^\e) COS0 + Bf^\e) sin 20. 



Taking the dot product of Eq. CIO with j yields: 



^ = 26f + (47r(sin^sin0 - g'')/)^'') + (C12) 

0(j) ^ ' 



1. Integrate Eq. \CT2\ 



fcC'^+i) = f\bf^'^ + (47r(sin^sin0 - g")/)^''^ + ^.fo^'^^rf^ + B^J;+^\e). 
Jo 

2. Impose periodicity of fei,^^^'' & integral constraint 

r 2bf^^^ + (47r(sin^sin0 - + £,6f+^)t^0 = 0, 



27r /»7r 



b^^^^Med(p = 0, 

'0 ^0 

to determine non-singular solutions for b'^^^\9). Note that calculation of 6^,^'' requires 
calculation of fc^^"* which requires calculation of f^^\ These calculations were performed 
using Maple, and files are available from the authors on request. 

We obtain the following expressions: 

hf = A(l -3cos2^)/108 
bW = B(l^\9) COS0 

6(2) = B^^'\9) + 5f )(^) sm0 + 5f )(^) cos 20 
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Taking the dot product of equation (CIO) with k yields: 



db 



47r/w cos ^ + Csbf\ 



(C13) 



which when combined with periodicity and the integral constrain yield the following expres- 
sions: 

1 



cos 9 



A sin(20) cos i 



= B^'^'iO) + B^''>{9) sin0 + B^''>{e) cos 20 



?(21)/ 



>(22)/ 



From equation (10) we can now compute the diffusion tensor correct to 0(l/cr^) : 

1 ,2 1 



J \2\ 

270 ' 



2430 
1 5 



1 , 2A2 41A^ 
(6 



243 25515' 



6 18a2 
810(7 



A^ 



(C14) 
(C15) 
(C16) 
(C17) 



all other entries are zero. When converting back to the global co-ordinates we note that 

i = e^, j = -e^, k = e^. (Cl^ 
and so with respect to er,e^, unit vectors the diffusion tensor is given by 



D 




















V 


a 





1 _ ^ 

6 ct2 



.d2 







ds + ^ J 



(C19) 



where 



di 



- H A^ 

3 270 

810 
A2 



2430 



dd = 6 



2A2 



41A^ 



243 25515 



-A^ 



(C20) 
(C21) 
(C22) 
(C23) 
(C24) 
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Appendix D: Fokker-Planck calculation of diffusion 

For the Fokker-Planck approximation,'^ the diffusion tensor non-dimensionalised on 
/dr is given by 

= Td, [{p- q)'f{p)dp, (Dl) 
Jp 

where r is a directional correlation time estimated from experimental data. Although the 
quantity r may vary with both A and the shear, a, for simplicity it is typically assumed to 
be independent of the shear. Asymptotic results for the diffusion tensor for weak sheai'^ 
(cr ^ 1) and strong shear^ (cr ^ 1) are available. With this choice of non-dimensionalisation 
and using the notation of this paper, the a ^ 1 result correct to 0(o"^) is given by 

D'; = Tdr^, D"-; = Td /^~^''^\ , D'p' = TdrK2, (D2) 

A A 

where K2 and J2 are specified functions of A'^.The a 3> 1 result correct to 0(l/o"^) is given 
by 

To compare the Fokker-Planck approximation to the generalized Taylor method we choose 
r so that the two alternative calculations for the horizontal component of the diffusion agree 
when the shear is zero. Specifically, when a = the generalized Taylor method yields 
Dq = and thus for the horizontal component of diffusion to agree we take rdr = 
For the specific gyrotactic bias A = 2.2, this yields rdr = 0.36. Taking this value of r also 
provides a value for the vertical component of diffusion. Dp = TdrK2 = 0.056, which is 
only a slight deviation from the generalized Taylor method, Dq = ^ = 0.050. Clearly, by 
this careful choice of r, the Fokker-Planck and generalized Taylor dispersion methods should 
agree when the shear is weak. However, as the shear increases we expect the two theories to 
give diverging predictions because, for example, in the FP approach Dp' approaches ^rdr, 
whereas in the GTD approach, Dq tends to zero at large shear. 

As for the generalized Taylor method, we fit simple functions to the curves of diffusion 



against a obtained with the Fokker-Planck methocP'^. The specific functions were given by 
D'^Xa) = P{a; aJT, b^^, Wi^) = -(yP{a- a^^ b^.^) D|f (a) = P((t; a|f , b^^). (D4) 



where the rational function P(o"; a, b) is defined by equation 19 and the choice of a coefficients 
is described in table HIl 
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TABLE II. In order to obtain the simplest functional fits whilst ensuring the asymptotic results 



are satisfied, the a coefficients are as specified. 





ao 


a2 


04 


^F 


Ji 

A2 
K2J1 

(Ki.h-J2).h 


JiX urr 1 Ji urr 
5i<'i"4,F 3KiX"2,F 

7JiX uzz 1 Ji uzz 
45Ki"4,F 1 3KiX"2,F 




Ji urr 
3KiX"4:,F 

J\ urr 
3ii'iA"4,F 





Appendix E: CoefRcients of functional fits 

The fit coefficients for A = 2.2 for the mean swimming and diffusion are given by: 





ao 


a2 


04 


b2 


&4 




2.05 X 10" 


-1 


1.86 X 10- 


-2 





1.74 X 10- 


-1 


1.27 X 10" 


-2 


a^ 


5.7 X 10- 


1 


3.66 X 10- 


-2 





1.75 X 10- 


-1 


1.25 X 10- 


-2 


„rr 


9.30 X 10- 


-2 


1.11 X 10" 


-4 





1.19 X 10- 


-1 


1.63 X 10" 


-4 




5.00 X 10^ 


-2 


1.11 X 10- 


-1 


3.71 X 10-^ 


1.01 X 10- 


-1 


1.86 X 10" 


-2 


^rz 
»G 


9.17 X 10- 


~2 


1.56 X 10" 


-4 





2.81 X 10- 


-1 


2.62 X 10" 


-2 


^rr 

ap 


9.30 X 10 


-2 


5.73 X 10" 


-4 


1.85 X 10-3 


4.96 X 10- 


-2 


1.54 X 10- 


-2 


nZZ 

ap 


5.60 X 10 


-2 


3.23 X 10- 


-2 


1.70 X 10-5 


2.70 X 10- 


-1 


1.42 X 10- 


-4 


^rz 

ap 


1.58 X 10- 


-2 










9.61 X 10- 


-2 


7.88 X 10- 


-2 
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